単位四元数\(q = \cos{\theta/2} \times \mathbf{1} + \sin{\theta/2} \times (x \mathbf{i} + y \mathbf{j} + z \mathbf{k});x^2+y^2+z^2=1\) を用いると、 3次元空間のベクトル\(v=(v_x,v_y,v_z)\)を軸\((x,y,z)\)の周りに角\(\theta\)で回転してできるベクトル\(v'=(v_x',v_y',v_z')\)は
\[ v_x' \mathbf{i} + v_y' \mathbf{j} + v_z' \mathbf{k} = q \times (v_x \mathbf{i} + v_y \mathbf{j} + v_z \mathbf{k}) \bar{q} \] ただし、\(\bar{q}=conj(q) = \cos{\theta/2} \times \mathbf{1} - \sin{\theta/2} \times (x \mathbf{i} + y \mathbf{j} + z \mathbf{k})\)
で得られる。
library(onion)
theta <- pi/3
xyz <- c(1,0,0)
v <- c(0,1,0)
q <- cos(theta/2) + sin(theta/2)*(xyz[1]*Hi + xyz[2]*Hj + xyz[3]*Hk)
v.q <- v[1]*Hi + v[2]*Hj + v[3]*Hk
v. <- q * v.q * Conj(q)
v.
## [1]
## Re 0.0000000
## i 0.0000000
## j 0.5000000
## k 0.8660254
Mod(v.)
## [1] 1
acos(sum(v*c(i(v.),j(v.),k(v.))))/pi # 内積を計算して、角度を求める、1/3 pi である
## [1] 0.3333333
library(onion)
tmp <- rnorm(4)
tmp <- tmp/sqrt(sum(tmp^2))
q <- tmp[1] + tmp[2]*Hi + tmp[3] * Hj + tmp[4] * Hk
v <- rnorm(3)
v.q <- v[1]*Hi + v[2] * Hj + v[3] * Hk
v. <- q * v.q * Conj(q)
A <- Re(q)
B <- i(q)
C <- j(q)
D <- k(q)
R <- matrix(c(A^2+B^2-C^2-D^2,2*(-A*D+B*C),2*(A*C+B*D),2*(A*D+B*C),A^2-B^2+C^2-D^2,2*(-A*B+C*D),2*(-A*C+B*D),2*(A*B+C*D),A^2-B^2-C^2+D^2),byrow=TRUE,3,3)
R
## [,1] [,2] [,3]
## [1,] 0.2917292 0.4464161 0.84593543
## [2,] -0.6065648 -0.5974892 0.52448620
## [3,] 0.7395764 -0.6661226 0.09647524
R %*% matrix(v,ncol=1)
## [,1]
## [1,] 2.7646513
## [2,] 0.3634669
## [3,] 0.7176218
v.
## [1]
## Re 0.0000000
## i 2.7646513
## j 0.3634669
## k 0.7176218
my.rot.mat <- function(V,theta){
V <- V/sqrt(sum(V^2))
A <- cos(theta/2)
B <- sin(theta/2) * V[1]
C <- sin(theta/2) * V[2]
D <- sin(theta/2) * V[3]
R <- matrix(c(A^2+B^2-C^2-D^2,2*(-A*D+B*C),2*(A*C+B*D),2*(A*D+B*C),A^2-B^2+C^2-D^2,2*(-A*B+C*D),2*(-A*C+B*D),2*(A*B+C*D),A^2-B^2-C^2+D^2),byrow=TRUE,3,3)
return(R)
}
my.function.tri.rot <- function(A,N,theta,zigzag=FALSE){
psi <- pi/3
if(zigzag){
psi <- - pi / 3
}
R <- my.rot.mat(A,theta)
N.new <- R %*% matrix(N,ncol=1)
#R2 <- my.rot.mat(N.new,pi/3)
R2 <- my.rot.mat(N.new,psi)
A.new <- R2 %*% matrix(A,ncol=1)
return(list(N=N.new,A=A.new))
}
my.function.tri.rot.series <- function(A=c(1,0,0),N=c(0,0,1),thetas,zigzag=FALSE){
As <- Ns <- matrix(0,length(thetas)+1,length(A))
As[1,] <- A
Ns[1,] <- N
for(i in 1:length(thetas)){
if(i %% 2 == 0){
tmp <- my.function.tri.rot(As[i,],Ns[i,],thetas[i],zigzag=FALSE)
}else{
tmp <- my.function.tri.rot(As[i,],Ns[i,],thetas[i],zigzag=zigzag)
}
As[i+1,] <- tmp$A
Ns[i+1,] <- tmp$N
}
return(list(As=As,Ns=Ns,k=length(thetas)))
}
A1 <- c(1,0,0)
N1 <- c(0,0,1)
k <- 6
thetas <- rep(0,k) # 回転角は0ばかり
ANs <- list()
ANs[[1]] <- list(A=A1,N=N1)
for(i in 1:k){
tmp <- my.function.tri.rot(ANs[[i]]$A,ANs[[i]]$N,thetas[i])
ANs[[i+1]] <- list(A=tmp$A,N=tmp$N)
}
# 確かに、辺ベクトル、法線ベクトルとが元に戻っている
ANs[[1]][[1]] - ANs[[k+1]][[1]]
## [,1]
## [1,] -4.440892e-16
## [2,] 1.276756e-15
## [3,] 0.000000e+00
ANs[[1]][[2]] - ANs[[k+1]][[2]]
## [,1]
## [1,] 0
## [2,] 0
## [3,] 0
その角が幾つなのか、計算機的に算出してみる。
角度を細かく変化させ、一周してきて、辺ベクトルと法線ベクトルとが一致する角度(の0より大きい)最小値を探索する。
まずは、2個。
これは、2つの正三角形の表裏貼り合わせ。 「正解」は\(\pi\)。
最初の辺と周回してきた辺との内積が1で、最初の法線ベクトルと周回してきた法線ベクトルとの内積が1であるような、角度が、閉じさせる角度。
3個は正四面体、4個はピラミッド、5個は正二十面体の1頂点を巡る5面、そして6個は平面。
7個以上になると、ある等角度で接続していくと「余って」しまうので、以下の例では、2周して閉じる場合が算出されている。
グラフ表示は、面の数ごとに2つ描いている。
左側のパネルは、横軸に回転角(単位\(\pi\))で、1周して戻ってきた、辺ベクトルと最初の辺ベクトルとの内積、最初と1周後の法線ベクトルの内積とを縦軸に取っている。
両方が揃って、内積1になる角度が、うまく閉じた場合に相当する。
うまく閉じる角度のうち、最小のものが、求める角度。
右側のパネルは、辺ベクトル内積と法線ベクトル内積との2次元変化の様子を表す。点(1,1)に到達すれば、それは、うまく閉じることを意味する。
正三角形の個数が増えると、うまく閉じる角度が複数通り現れる様子が、点(1,1)に繰り返し曲線が戻ることから読み取れる。
my.tri.cycle <- function(k, s=3,numtheta = 1003,max.theta = 2 * pi * 1,sign.alt=FALSE){
theta.list <- seq(from=0,to=max.theta,length=numtheta)
cosAs <- cosNs <- rep(0,length(theta.list))
cosANs.sq <- cosAs
for(ii in 1:length(theta.list)){
thetas <- rep(theta.list[ii],k)
if(sign.alt){
thetas <- thetas * (-1)^(1:length(thetas)) * (-1)
}
ANs <- list()
ANs[[1]] <- list(A=A1,N=N1)
for(i in 1:k){
tmp <- my.function.tri.rot(ANs[[i]]$A,ANs[[i]]$N,thetas[i])
ANs[[i+1]] <- list(A=tmp$A,N=tmp$N)
}
cosAs[ii] <- sum(ANs[[1]][[1]] * ANs[[k+1]][[1]])
cosNs[ii] <- sum(ANs[[1]][[2]] * ANs[[k+1]][[2]])
}
#theta.list[which((abs(cosAs-1) + abs(cosNs-1) < 10^(-4)) & ((abs(cosAs-1) + abs(cosNs-1) == min(abs(cosAs-1) + abs(cosNs-1)))))]
cosANs <- cbind(cosAs,cosNs)
cosANs.sum <- apply(cosANs,1,sum)
best.theta <- theta.list[which(cosANs.sum > 2-10^(-s))][1]
par(mfcol=c(1,2))
matplot(theta.list/pi,cbind(cosAs,cosNs),type="l",xlab="angle(unit=pi)",ylab="cosines of edge vectors and normal vectors",main = paste("k=",k))
abline(h=1,col=3)
abline(v=best.theta/pi,col=3)
plot(cosAs,cosNs,type="l",xlab="cos(edge_vecs)",ylab="cos(norm_vecs)")
points(c(1,1),pch=20,cex=3,col=2)
par(mfcol=c(1,1))
return(list(theta=best.theta,cosAs=cosAs,cosNs=cosNs,theta.list=theta.list,sign.alt=sign.alt))
}
ks <- 2:30
outs <- list()
for(i in 1:length(ks)){
outs[[i]] <- my.tri.cycle(k=ks[i])
}
三角形の枚数と均等角との関係は次のようになる。
estimated.theta <- sapply(outs, function(x){return(x$theta)})/pi
plot(ks,estimated.theta,ylab="angle, unit=PI")
#abline(h=c(1,1/2,1/3,1/4))
正三角形の枚数が2、4、8枚の場合を重ねてプロットするとわかるように、枚数の約数の「周回角度」が含まれる。
matplot(cbind(outs[[1]]$cosAs,outs[[3]]$cosAs,outs[[7]]$cosAs),type="l")
2枚、3枚、6枚では?
matplot(cbind(outs[[1]]$cosAs,outs[[2]]$cosAs,outs[[5]]$cosAs),type="l")
5枚、6枚、10枚、30枚では?
matplot(cbind(outs[[4]]$cosAs,outs[[5]]$cosAs,outs[[9]]$cosAs,outs[[29]]$cosAs),type="l")
この「角度」と枚数とその約数との関係に、ある種の分子分母関係があることがわかる。
上の例では、7枚以上においても、均等角を算出した。
しかしながら、この角では、三角形同士が交叉することとなり、3次元空間の閉多面体としては不適当である。
したがって、7枚以上の場合には、等角度ではあるが、角度の正負が交互になるような均質性を仮定して、閉じさせることを考える。
実際に、実施してみると、奇数枚の場合には、正負が反転しない場合が発生し、そのような閉じ方が成立しないことがわかる。
ks2 <- 7:30
outs.alt <- list()
for(i in 1:length(ks2)){
outs.alt[[i]] <- my.tri.cycle(k=ks2[i],sign.alt=TRUE) # sign.alt=TRUEにより、角の正負を交互にする
}
算出された角を、上段に偶数枚の場合、下段に奇数枚の場合を書くと、偶数枚の場合には値がありるが、奇数枚の場合には、値が得られないことがわかる。
outs.alt.ang <- sapply(outs.alt,function(x){return(x$theta)})
matrix(outs.alt.ang,nrow=2)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] NA NA NA NA NA NA NA NA NA
## [2,] 1.385812 1.805945 0 1.034656 1.392083 0 0.8590782 1.178881 0
## [,10] [,11] [,12]
## [1,] NA NA NA
## [2,] 0.7524773 1.040927 0
plot(matrix(ks2,nrow=2)[2,],matrix(outs.alt.ang,nrow=2)[2,])
角度の符号反転の場合には、枚数の約数関係に、反転無しの場合のような規則は存在しないようである。
以下に、8枚、16枚の関係を示す。
matplot(cbind(outs.alt[[2]]$cosAs,outs.alt[[10]]$cosAs),type="l")
正三角形を1点を中心にk枚貼り合わせて、元に戻るようにする。
ただし、k = 6t + s の時には、s = 0 ならば、\(2 \times t \times \pi\)の角をk等分することとし、s = 1,2,…,5ならば、\(2 \times (t+1) \times \pi\) の角をk等分することとする。
計算してみる。
\[ \phi = \frac{2 t \pi}{k}\\ L = \frac{1}{2 \sin{\frac{\phi}{2}}\\ \cos{\theta} = \frac{4}{3} * (\frac{1}{2 L^2} - \frac{5}{4}) \]
関数にする。
# k は三角形の枚数
my.tri.angle <- function(k){
# k = t * 6 の時は、角度0
# k = t * 6 + s (s<6) の時は、2 pi * (t+1) の角をk等分する
six.info <- (k-1)%/%6
# six.info + 1 周させる
# 1つの三角形が稼ぐべき、二次元平面的な角度
phi <- 2 * pi * (six.info + 1) / k
# 中心頂点から見下ろした時に、正三角形の周辺頂点の原点からの距離(2D平面上の距離)
L <- 1/(2*sin(phi/2))
# 幾何的な関係を色々やると、正三角形の周辺頂点を一つ飛ばしにした2頂点と
# その間にある周辺頂点と中心頂点の中点とをそれぞれ結び、その2つの線分ベクトルのcos(theta)を計算すると以下の式となる
costheta <- 4/3 * (1/(2*L^2) - 5/4)
# 計算誤差のためにちょっと工夫をして
if(costheta > 1) costheta <- 1
if(costheta < -1) costheta <- -1
theta <- acos(costheta)
# 最終的に、2つの三角形のなす角を、pi-thetaとして返す
return(pi - theta)
}
多数の角を評価して推定した結果と一致することを確かめる。
# ks <- 2:30
t <- rep(0,length(ks))
for(i in 1:length(ks)){
t[i] <- my.tri.angle(ks[i])
}
plot(t)
plot(estimated.theta, t/pi)
abline(0,1,col=2)
\(k = 2p\)枚の正三角形を用いる。
外周頂点のz座標は絶対値は等しく、正負が交互になるから、ある2枚の正三角形は、原点(0,0,0)を中心として、
(x,-y,-z), (x’,0,z),(x,y,-z)の3点となる。
こららが満足する条件は \[ x^2 + y^2 + z^2 = x^2 + (-y)^2 + (-z)^2 = x^2 + y^2 + (-z)^2 = 1\\ x'^2 + z^2 = 1\\ (x-x')^2 + (-y)^2 + (2z)^2 = 1\\ \] また、外周辺の中点のz座標は0であり、隣り合う、中心と外周辺中点を結ぶベクトルとのなす角は、\(\frac{2\pi}{k}\)である。
中点2つは\((\frac{x+x'}{2},-\frac{y}{2},0),(\frac{x+x'}{2},\frac{y}{2},0)\)であるから \[ \frac{(\frac{x+x'}{2})^2 - (\frac{y}{2})^2}{(\frac{x+x'}{2})^2 + (\frac{y}{2})^2} = \cos{\frac{2\pi}{k}} \] これらを解くと \[ x' = \sqrt{\frac{3}{2}\frac{1}{1+\cos{\frac{2\pi}{k}}}}\\ z = \sqrt{1-x'^2}\\ x = x' \cos{\frac{2\pi}{k}}\\ y = x' \sqrt{(1+\cos{\frac{2\pi}{k}})(1-\cos{\frac{2\pi}{k}})} \] 2枚の三角形がなす角\(\theta\)は (x-x’/2,-y,-3/2z) と(x-x’/2,y,-3/2z)とがなす角の外角であるから、
\[ \cos{(\pi-\theta)} = \frac{(x-\frac{x'}{2})^2-y^2+9/4z^2}{(x-\frac{x'}{2})^2+y^2+9/4z^2} \]
これを関数にして
my.updown.ang <- function(k){
C <- cos(2*pi/k)
x. <- sqrt(3/(2 * (1+C)))
z <- sqrt(1-x.^2)
x <- x. * C
y <- x. * sqrt((1+C)*(1-C))
tmp.cos <- ((x-x./2)^2 - y^2 + 9/4*z^2)/((x-x./2)^2 + y^2 + 9/4*z^2)
theta <- pi - acos(tmp.cos)
return(theta)
}
out.alt.2 <- rep(0,length(ks2))
for(i in 1:length(ks2)){
out.alt.2[i] <- my.updown.ang(ks2[i])
}
plot(ks2,out.alt.2)
角度をパラメタ化して探索した場合は、2周以上しても良いので、なるべく小さい角で元に戻る値を求めているため、幾何的算出結果と合致しない。
outs.alt.ang
## [1] NA 1.3858123 NA 1.8059455 NA 0.0000000 NA
## [8] 1.0346563 NA 1.3920830 NA 0.0000000 NA 0.8590782
## [15] NA 1.1788811 NA 0.0000000 NA 0.7524773 NA
## [22] 1.0409269 NA 0.0000000
out.alt.2
## [1] 1.040208 1.398370 1.635018 1.809114 1.944444 2.053416 2.143402 2.219151
## [9] 2.283897 2.339937 2.388956 2.432220 2.470703 2.505168 2.536222 2.564352
## [17] 2.589958 2.613368 2.634855 2.654649 2.672944 2.689905 2.705673 2.720372
枚数が偶数の場合について、悉皆探索最小角と、幾何的産出角の関係を示す。
plot(matrix(outs.alt.ang,nrow=2)[2,],matrix(out.alt.2,nrow=2)[2,])
合致関係は、以下のようなグラフとして表示するとわかりやすい。
また、奇数枚の場合がなぜ不適当かも、視覚的に表現されている。
for(i in 1:length(ks2)){
tmp <- cbind(outs.alt[[i]]$cosAs,outs.alt[[i]]$cosNs)
matplot(outs.alt[[i]]$theta.list/pi,tmp,type="l")
abline(v = out.alt.2[i]/pi,col=4)
}
今、k枚の正三角形が、 \[ (\theta_1,\theta_2,...,\theta_k) \] で閉じているとする。
k-1枚目とk枚目との間に2枚を\(\pi\)でピタリと押しつぶした状態で割り込ませ、k+2枚にすることができる。 その時にできる、長さk+2の角度列は
\[ (\theta_1,...,\theta_{k-1},-\pi + \phi, \pi, \theta_k - \phi) \]
5枚で均等曲げを実行してみる。
k <- 5
theta.k <- t[k-1]
thetas <- rep(theta.k,k)
outk <- my.function.tri.rot.series(thetas=thetas)
range(outk$As[1,] -outk$As[k+1,])
## [1] -2.867186e-16 -1.110223e-16
range(outk$Ns[1,] -outk$Ns[k+1,])
## [1] 0.000000e+00 2.220446e-16
my.insert.two.tris <- function(theta,phi = 0.1){
k <- length(theta)
ret <- c(theta[1:(k-1)],-pi + phi,pi,theta[k]-phi)
return(ret)
}
thetas2 <- my.insert.two.tris(thetas)
k2 <- length(thetas2)
outk2 <- my.function.tri.rot.series(thetas=thetas2)
range(outk2$As[1,] -outk2$As[k2+1,])
## [1] -4.86452e-16 0.00000e+00
range(outk2$Ns[1,] -outk2$Ns[k2+1,])
## [1] -6.661338e-16 3.330669e-16
2枚貼り合わせ三角形を挿入して角度ベクトルを不均一にした。それを改変して、閉じる・閉じない条件がどうなるか・・・
角度の順序を入れ替えると、閉じなくなる。。。。単純ではない。そりゃそうか。
thetas2.shuffle <- sample(thetas2)
outk3 <- my.function.tri.rot.series(thetas=thetas2.shuffle)
range(outk3$As[1,] -outk3$As[k2+1,])
## [1] -5.178593e-16 -2.220446e-16
range(outk3$Ns[1,] -outk3$Ns[k2+1,])
## [1] 0.000000e+00 3.816392e-16
1頂点周りにk枚の正三角形が並んで閉じる例を検討してきた。
今度は、正三角形がジグザグ路・擬直線周回路をなす場合について検討する。
すでに作成した関数の一部改変で可能である。
4枚で一周する時は、正四面体様の形になるので、3枚を1頂点周りにとった場合と同じ角度で閉じることに注意する。
6枚の場合に\(\pi\)でも辺ベクトルと法線ベクトルとが一致するが、ベクトルの向きとしてそのようになるだけであって、閉じるわけではないことに注意する。
my.tri.zigzag.cycle <- function(k, s=3,numtheta = 1003,max.theta = 2 * pi * 1,sign.alt=FALSE){
theta.list <- seq(from=0,to=max.theta,length=numtheta)
cosAs <- cosNs <- rep(0,length(theta.list))
cosANs.sq <- cosAs
for(ii in 1:length(theta.list)){
thetas <- rep(theta.list[ii],k)
if(sign.alt){
thetas <- thetas * (-1)^(1:length(thetas)) * (-1)
}
ANs <- list()
ANs[[1]] <- list(A=A1,N=N1)
for(i in 1:k){
if(i %% 2 == 0){
tmp <- my.function.tri.rot(ANs[[i]]$A,ANs[[i]]$N,thetas[i],zigzag=FALSE)
}else{
tmp <- my.function.tri.rot(ANs[[i]]$A,ANs[[i]]$N,thetas[i],zigzag=TRUE)
}
ANs[[i+1]] <- list(A=tmp$A,N=tmp$N)
}
cosAs[ii] <- sum(ANs[[1]][[1]] * ANs[[k+1]][[1]])
cosNs[ii] <- sum(ANs[[1]][[2]] * ANs[[k+1]][[2]])
}
#theta.list[which((abs(cosAs-1) + abs(cosNs-1) < 10^(-4)) & ((abs(cosAs-1) + abs(cosNs-1) == min(abs(cosAs-1) + abs(cosNs-1)))))]
cosANs <- cbind(cosAs,cosNs)
cosANs.sum <- apply(cosANs,1,sum)
best.theta <- theta.list[which(cosANs.sum > 2-10^(-s))][1]
par(mfcol=c(1,2))
matplot(theta.list/pi,cbind(cosAs,cosNs),type="l",xlab="angle(unit=pi)",ylab="cosines of edge vectors and normal vectors",main = paste("k=",k))
abline(h=1,col=3)
abline(v=best.theta/pi,col=3)
plot(cosAs,cosNs,type="l",xlab="cos(edge_vecs)",ylab="cos(norm_vecs)")
points(c(1,1),pch=20,cex=3,col=2)
par(mfcol=c(1,1))
return(list(theta=best.theta,cosAs=cosAs,cosNs=cosNs,theta.list=theta.list,sign.alt=sign.alt))
}
ks <- 2:30
outs.zigzag <- list()
for(i in 1:length(ks)){
outs.zigzag[[i]] <- my.tri.zigzag.cycle(k=ks[i])
}
4枚と6枚と12枚では、12枚の方が4枚の角度で三巡する場合、6枚の角度で二巡する場合があることが、以下のプロットからわかる
matplot(cbind(outs.zigzag[[3]]$cosAs,outs.zigzag[[5]]$cosAs,outs.zigzag[[11]]$cosAs),type="l")
2個の正三角形の頂点座標は次のように表される。
\(X = (L\cos{\frac{2\pi}{k}},-L\sin{\frac{2\pi}{k}},0)\), \(Y = (L\cos{\frac{2\pi}{k}},L\sin{\frac{2\pi}{k}},0)\), \(Z = (L,0,z)\), \(W = (L\cos{\frac{4\pi}{k}},L\sin{\frac{4\pi}{k}},z)\)
XYZ,YZWが正三角形である。
2枚の正三角形のなす角は、X–(Y,Zの中点)、(Y,Zの中点)–W の2ベクトルのなす角である。
正三角形である条件から、XY,YZ,XZ,YW,ZWの長さが1である。 したがって
\[ 2L\sin{\frac{2\pi}{k}} = 1\\ (L\cos{\frac{2\pi}{k}} - L)^2 + (L\sin{\frac{2\pi}{k}})^2 + z^ 2=1 \]
これを解いて \[ L = \frac{1}{2\sin{\frac{2\pi}{k}}}\\ z = \sqrt{1-\frac{1-\cos{\frac{2\pi}{k}}}{2\sin{\frac{2\pi}{k}}^2}} \] 関数を作る。
my.zigzag.ang <- function(k){
ang <- 2*pi/k
L <- 1/(2*sin(ang))
z <- sqrt(1-(1-cos(ang))/(2*sin(ang)^2))
X <- c(L*cos(ang),-L*sin(ang),0)
Y <- c(L*cos(ang),L*sin(ang),0)
Z <- c(L,0,z)
W <- c(L*cos(ang*2),L*sin(ang*2),z)
YZ <- (Y + Z)/2
X.YZ <- X - YZ
W.YZ <- W - YZ
tmpcos <- sum(X.YZ*W.YZ)/(sqrt(sum(X.YZ^2))*sqrt(sum(W.YZ^2)))
theta <- pi - acos(tmpcos)
return(theta)
}
計算しておく。
ang.zigzag2 <- rep(0,length(ks))
for(i in 1:length(ks)){
ang.zigzag2[i] <- my.zigzag.ang(ks[i])
}
## Warning in sqrt(1 - (1 - cos(ang))/(2 * sin(ang)^2)): 計算結果が NaN になりまし
## た
plot(ks,ang.zigzag2)
悉皆探索の結果との一致を確認する。kが偶数の時にきちんと計算が合っていることがわかる
for(i in 1:length(ks)){
tmp <- cbind(outs.zigzag[[i]]$cosAs,outs.zigzag[[i]]$cosNs)
matplot(outs.zigzag[[i]]$theta.list/pi,tmp,type="l",main=ks[i])
abline(v = ang.zigzag2[i]/pi,col=4)
}
## まとめ
ひとかたまりの何かを整数で等分したものを、単位分数と呼ぶことにする。
\(\frac{1}{k}; k = 1,2,...\) が単位分数。
正三角形を3次元空間に並べることによって、複数の単位分数的角度というものが定まる。
このような、複数種類ある正三角形の3次元空間配置に伴う単位分数的角度のみを使って、正三角形メッシュはできているのではないか、という仮説を立てた。
これは正しいか?